home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / c_correlate.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  129 lines
  1. ; $Id: c_correlate.pro,v 1.8 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1995-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       C_CORRELATE
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the cross correlation Pxy(L) or cross 
  11. ;       covariance Rxy(L) of two sample populations X and Y as a function
  12. ;       of the lag (L).
  13. ;
  14. ; CATEGORY:
  15. ;       Statistics.
  16. ;
  17. ; CALLING SEQUENCE:
  18. ;       Result = C_correlate(X, Y, Lag)
  19. ;
  20. ; INPUTS:
  21. ;       X:    An n-element vector of type integer, float or double.
  22. ;
  23. ;       Y:    An n-element vector of type integer, float or double.
  24. ;
  25. ;     LAG:    A scalar or n-element vector, in the interval [-(n-2), (n-2)],
  26. ;             of type integer that specifies the absolute distance(s) between
  27. ;             indexed elements of X.
  28. ;
  29. ; KEYWORD PARAMETERS:
  30. ;       COVARIANCE:    If set to a non-zero value, the sample cross 
  31. ;                      covariance is computed.
  32. ;
  33. ;       DOUBLE:        If set to a non-zero value, computations are done in
  34. ;                      double precision arithmetic.
  35. ;
  36. ; EXAMPLE
  37. ;       Define two n-element sample populations.
  38. ;         x = [3.73, 3.67, 3.77, 3.83, 4.67, 5.87, 6.70, 6.97, 6.40, 5.57]
  39. ;         y = [2.31, 2.76, 3.02, 3.13, 3.72, 3.88, 3.97, 4.39, 4.34, 3.95]
  40. ;
  41. ;       Compute the cross correlation of X and Y for LAG = -5, 0, 1, 5, 6, 7
  42. ;         lag = [-5, 0, 1, 5, 6, 7]
  43. ;         result = c_correlate(x, y, lag)
  44. ;
  45. ;       The result should be:
  46. ;         [-0.428246, 0.914755, 0.674547, -0.405140, -0.403100, -0.339685]
  47. ;
  48. ; PROCEDURE:
  49. ;       See computational formula published in IDL manual.
  50. ;
  51. ; REFERENCE:
  52. ;       INTRODUCTION TO STATISTICAL TIME SERIES
  53. ;       Wayne A. Fuller
  54. ;       ISBN 0-471-28715-6
  55. ;
  56. ; MODIFICATION HISTORY:
  57. ;       Written by:  GGS, RSI, October 1994
  58. ;       Modified:    GGS, RSI, August 1995
  59. ;                    Corrected a condition which excluded the last term of the
  60. ;                    time-series.
  61. ;       Modified:    GGS, RSI, April 1996
  62. ;                    Simplified CROSS_COV function. Added DOUBLE keyword.
  63. ;                    Modified keyword checking and use of double precision.
  64. ;-
  65.  
  66. FUNCTION Cross_Cov, X, Y, M, nX, Double = Double
  67.   ;Sample cross covariance function.
  68.  
  69.   Xmean = TOTAL(X, Double = Double) / nX
  70.   Ymean = TOTAL(Y, Double = Double) / nX
  71.   RETURN, TOTAL((X[0:nX - M - 1L] - Xmean) * (Y[M:nX - 1L] - Ymean), $
  72.                  Double = Double)
  73.  
  74. END
  75.  
  76. FUNCTION C_Correlate, X, Y, Lag, Covariance = Covariance, Double = Double
  77.  
  78.   ;Compute the sample cross correlation or cross covariance of 
  79.   ;(Xt, Xt+l) and (Yt, Yt+l) as a function of the lag (l).
  80.  
  81.   ON_ERROR, 2
  82.  
  83.   TypeX = SIZE(X)
  84.   TypeY = SIZE(Y)
  85.   nX = TypeX[TypeX[0]+2]
  86.   nY = TypeY[TypeY[0]+2]
  87.  
  88.   if nX ne nY then $
  89.     MESSAGE, "X and Y arrays must have the same number of elements."
  90.  
  91.   ;Check length.
  92.   if nX lt 2 then $
  93.     MESSAGE, "X and Y arrays must contain 2 or more elements."
  94.  
  95.   ;If the DOUBLE keyword is not set then the internal precision and
  96.   ;result are identical to the type of input.
  97.   if N_ELEMENTS(Double) eq 0 then $
  98.     Double = (TypeX[TypeX[0]+1] eq 5 or TypeY[TypeY[0]+1] eq 5)
  99.  
  100.   nLag = N_ELEMENTS(Lag)
  101.  
  102.   if nLag eq 1 then Lag = [Lag] ;Create a 1-element vector.
  103.  
  104.   if Double eq 0 then Cross = FLTARR(nLag) else Cross = DBLARR(nLag)
  105.  
  106.   if KEYWORD_SET(Covariance) eq 0 then begin ;Compute Cross Correlation.
  107.     for k = 0, nLag-1 do begin
  108.       if Lag[k] ge 0 then $
  109.         Cross[k] = Cross_Cov(x, y, Lag[k], nX, Double = Double) / $
  110.                    SQRT(Cross_Cov(X, X, 0L, nX, Double = Double) * $
  111.                    Cross_Cov(Y, Y, 0L, nY, Double = Double)) $
  112.       else Cross[k] = Cross_Cov(Y, X, ABS(Lag[k]), nY, Double = Double) / $
  113.                    SQRT(Cross_Cov(X, X, 0L, nX, Double = Double) * $
  114.                    Cross_Cov(Y, Y, 0L, nY, Double = Double))
  115.     endfor
  116.   endif else begin ;Compute Cross Covariance.
  117.     for k = 0, nLag-1 do begin
  118.       if Lag[k] ge 0 then $
  119.         Cross[k] = Cross_Cov(X, Y, Lag[k], nX, Double = Double) / nX $
  120.       else Cross[k] = Cross_Cov(Y, X, ABS(Lag[k]), nX, Double = Double) / nX
  121.     endfor
  122.   endelse
  123.  
  124.   if Double eq 0 then RETURN, FLOAT(Cross) else $
  125.   RETURN, Cross
  126.  
  127. END
  128.   
  129.